Black hole formation through fragmentation of toroidal polytropes 
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We investigate new paths to black hole formation by considering the general relativistic evolution of a dif- 
ferentially rotating polytrope with toroidal shape. We find that this polytrope is unstable to nonaxisymmetric 
modes, which leads to a fragmentation into self-gravitating, collapsing components. In the case of one such 
fragment, we apply a simplified adaptive mesh refinement technique to follow the evolution to the formation of 
an apparent horizon centered on the fragment. This is the first study of the one-armed instability in full general 
relativity. 



The formation of black holes from neutron stars, iron cores 
or supermassive stars is expected to be associated with a char- 
acteristic gravitational wave signal which may give informa- 
tion about the collapse dynamics and the physical environ- 
ment of such objects. Therefore, and given that gravitational 
wave detectors are already taking data or are coming online, 
it is of prime importance to understand the dynamical features 
of the gravitational collapse of hydrodynamical systems. 

The prototypical model of stellar collapse is an equilibrium 
polytrope subject to a radial or quasi-radial perturbation grow- 
ing on a dynamical timescale. In spherical symmetry, every 
general relativistic polytrope with index N = 3 is unstable to 
radial oscillations 1 1 ] - in turn, there exists a critical N c < 3 
for which the star is marginally stable. Without spherical sym- 
metry, rotation can increase this critical value again |2[. The 
black hole formation from the collapse of uniformly and dif- 
ferentially rotating polytropes induced by this instability is a 
well-investigated phenomenon, either with restriction to ax- 
isymmetry S 011 II 0111] or without O [H E d . In 
the gauge choices usually employed, the dynamical behaviour 
of the system shows a radial contraction of the star, accompa- 
nied by the formation of an apparent horizon at late times. 

Black hole formation from a dynamically unstable nonax- 
isymmetric mode, however, has not been modelled so far. Sce- 
narios range from the development of a bar mode, subsequent 
transport of angular momentum into the shell and collapse 
of the central object, to fragmentation and off-center produc- 
tion of one or several black holes. In Newtonian theory, in- 
stabilities and fragmentation have received considerable at- 
tention, specifically in the context of binary formation from 
protostellar disks (e.g. 

giee Gam 

and references 

therein) and compact object production in stellar core collapse 
(e.g. 1191 l20l 12. ill and references therein). In II ill , the authors 
also report signatures of an m = 4 fragmentation behaviour 
in the collapse and centrifugal bounce of an N = 1 polytrope, 
but could not determine the final state due to resolution issues. 

The cooling evolution of supermassive stars can be approx- 
imately described by the N = 3 mass-shedding sequence 
when the angular momentum transport timescales are short 
compared to the cooling timescale 1 22], so that uniform ro- 
tation is enforced. This sequence has a turning point for the 



onset of a quasi-radial instability, and numerical experiments 
confirm that the collapse remains axisymmetric 1 23 ] . If the 
star is differentially rotating, the cooling sequence is less con- 
strained and mig ht end in a transition to nonaxisymmetric in- 
stability 124 125L 26]. The canonical expectation that a super- 
massive star produces one central black hole with a low-mass 
accretion disk might thus not be appropriate for differentially 
rotating configurations. 

In this letter, we consider the production of a black hole 
through the fragmentation of a general relativistic polytrope. 
We focus on N — 3 polytropes, which are associated with pre- 
collapse cores of massive stars or supermassive stars. The soft 
equation of state also enhances the instability of the fragments 
compared to the common choice N = 1 for neutron stars. To 
represent this process accurately on a grid, we make use of an 
adaptive mesh refinement technique, since a possibly highly 
deformed apparent horizon needs to be located in some region 
of the domain which is unknown in advance. 

The recent investigation of the collapse of differentially ro- 
tating supermassive stars by Saijo 1 27] was based on a se- 
quence of relativistic N ~ 3 polytropes with a parameterized 
rotation law of the commonly used form j($Y) = A 2 (fi c — fi), 
where £! c is the angular velocity at the center, and the param- 
eter A specifies the degree of differential rotation (A — > oo 
is uniform rotation). The sequence selected was constrained 
by a constant central density p c = 3.38 x 1CP 6 in units 
K = G = c = 1, and the choice A/r c = 1/3, where r c 
denotes the equatorial coordinate radius. 

To examine the indirect collapse by fragmentation of a 
polytrope with toroidal shape, we choose a model with the 
same central density as in Saijo's 12711 models, but with a ratio 
of polar to equatorial coordinate radius r p /r = 0.24. The ra- 
tio of rotational kinetic energy to gravitational binding energy 
is T/\W\ = 0.227. While the critical limit for the dynami- 
cal f-mode instability in uniform density, uniformly rotating 
Maclaurin spheroids is (T/|W|) dyn = 0.2738 (e.g. jH), re- 
cent investigations of the stability of soft (N ~ 3) differen- 
tially rotating polytropes in Newtonian gravity 1 29. l30l.l3 lLl32ll 
have shown that the Maclaurin approximation is inappropriate 
for such systems, and generally find the critical (r/|W|) d 
to be below the Maclaurin value. Consequently, the toroid- 
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like star considered in this study might be unstable to nonax- 
isymmetric perturbations. Here we present the first investiga- 
tion of this instability in full general relativity, showing that 
relativistic effects are significant for the final outcome, as we 
observe that black holes can be produced. 

All simulations have been performed in full general rela- 
tivity. The only assumption on symmetry is a reflection in- 
variance with respect to the equatorial plane of the star. The 
gauge freedom is fixed by the generalized 1+log slicing con- 
dition for the lapse function 13311 with f(a) =2/ a, and by the 
hyperbolic-type condition suggested in ||8|] for the shift vector. 

The computational framework is the Cactus code 
( |www . cactuscode . org) , which also provides a module 
to solve the geometric part of the field equations in the 
well-known BSSN form |H Q. In addition, the 
Carpet driver 113 711 is used for mesh refinement in Cactus. 
The hydrodynamics part of the field equations is evolved 
using the high-resolution shock-capturing PPM-Marquina 
implementation in the Whisky module 111 2113811 . and a gamma 
law equation of state (P = pe/N). We are thus using a 
set of well-tested tools to evolve the general relativistic 
hydrodynamics equations. 

To numerically construct the axisymmetric initial model de- 
scribed in the introduction, we use the RNS initial data solver 
ll39ll with a radial resolution of 601 and an angular resolution 
of 301 points. With the parameters described above, the model 
has toroid-like structure, with an off-center density maximum, 
but a non-zero central density. After mapping the model to the 
hierarchy of Cartesian grids provided by Carpet, a small per- 
turbation of the form 



p(x) ->■ p(x) 
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is applied with A m = 0, 1 and A = £\ Ai. In addition, the 
polytropic constant K is reduced by 0. 1% to induce collapse if 
the model is radially unstable. After perturbing the model, the 
constraint equations are not solved again, since the amplitude 
B is chosen such that the violation of the constraints by the 
initial perturbation is about an order of magnitude smaller than 
that caused by the systematic error induced by the m = 4 
symmetry of the Cartesian grid. 

For most simulations, a fixed box-in-box mesh refinement 
with 5 levels is used to accurately resolve the central high- 
density ring. The three innermost grids cover the star, while 
the two outermost ones push the outer boundaries to 6.4r e . 
The typical resolution used was 65 x 65 x 33 per grid patch, 
leading to a central resolution of 5 X rj 10~ 2 r c . However, 
runs with 89 x 89 x 45, 97 x 97 x 52 and 131 x 131 x 65 
points per grid patch were also performed to test the resolution 
dependence of the solution; for the last setup, a simulation 
with a uniform grid setup would need to cover the equatorial 
plane of the star alone with 320 grid points to achieve the same 
central resolution. 

To determine the amplitude of a specific mode in the equa- 
torial plane, we perform a projection onto Fourier modes at 
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FIG. 1 : Time evolution of the mode amplitudes in the standard grid 
setup N. The amplitudes are obtained from a Fourier decomposition 
of the density profile on the equatorial plane circle at w = 0.25r c , 
the initial radius of highest density. Shown are the m = 1 (thick solid 
line), m — 2 (thick dashed line), m — 3 (thin dash-dotted line) and 
m = 4 (thin dotted line) mode amplitudes. The time is normalized 
to the dynamical timescale tu = r c r c /M. 



certain coordinate radii I40I EUl . Care must be taken in in- 
terpreting the results as soon as the system starts to deviate 
significantly from axisymmetry, since the interpretation of the 
projection curve as a circle assumes to be a Killing vector. 
In addition to the Fourier extraction, we monitor the evolution 
of the rest mass and the ADM constraints. 

In Table |l] we have listed the parameters for the different 
simulation runs. For model N, the evolution of the moduli of 
the equatorial Fourier components at the initial radius of high- 
est density is shown in Fig ^ It is evident that, initially, the 
m = 4 component (thin dotted line) induced by our Carte- 
sian grid is dominant. However, the star is unstable to m = 1 
(thick solid line) and m — 2 (thick dashed line), and these 
modes consequently grow into the nonlinear regime, their e- 
folding times being rather close. 

To test the dependence of the results on the grid resolution, 
we have evolved the same initial data with different grid pa- 
rameters as listed in Table U Runs HI, H2 and H3 are high 
resolution versions of N. The setups II and 12 test the depen- 
dence on the amplitude B of the initial data perturbation, with 

11 using an amplitude 10 times smaller than in setup N, and 

12 using B = 0. Ml and M2 are variants of N, where only a 
specific unstable mode is imposed. 

The rest mass is conserved numerically within 1.8% for 
setup N and to within 0.2% for setup H3. An approximate 
measurement of the e-folding times and mode frequencies can 
be obtained within an error of 5 — 10% related to ambiguities 



TABLE I: Parameters for the simulation runs. Different setups were 
chosen to confirm the results, including resolution tests, changes in 
the initial perturbation, and two setups with adaptive mesh refine- 
ment (AMR) to investigate black hole formation. 





Patch 


Refinement 


Perturbation 




Setup 


resolution 


levels 


Modes 


Amplitude 


AMR 


N 


65 x 65 x 33 


5 


m = 1 — 4 


B = 10" 3 


no 


HI 


89 x 89 x 45 


5 


m = 1 — 4 


B = 10" 3 


no 


H2 


97 x 97 x 52 


5 


m = 1 — 4 


B = 10~ 3 


no 


H3 


131 x 131 x 65 


5 


m = 1 — 4 


B = 1(T 3 


no 


11 


65 x 65 x 33 


5 


m = 1 — 4 


B = 1(T 4 


no 


12 


65 x 65 x 33 


5 


none 


5 = 


no 


Ml 


65 x 65 x 33 


5-12 


m — 1 


B = 10" 3 


yes 


M2 


65 x 65 x 33 


5-12 


m = 2 


B = 10" 3 


yes 



in defining the interval of extraction. All setups show consis- 
tent results within this uncertainty. In units of the dynamical 
timescale, which is defined here as £d = r c y/r c /M, the e- 
folding times are w 0.93£d for m = 1, and w 0.84£d for 
m = 2, respectively. Mode frequencies are » 3.05/£d for 
m = 1 and w 3.31/£d for m = 2, respectively. 

To establish whether a black hole is formed by a fragment 
it is necessary to cover the fragment with significantly more 
resolution than affordable by fixed mesh refinement. Hence 
we have implemented a simplified adaptive mesh refinement 
scheme to follow the system to black hole formation: In this 
scheme, a tracking function, here provided by the location of 
a density maximum, is used to construct a locally fixed hierar- 
chy of grids moving with the fragment. Additional refinement 
levels are switched on during contraction, until an apparent 
horizon is found. 

Since the e-folding times for m = 1 and m = 2 turn out 
to be close, the number and interaction behaviour of the frag- 
ments in the non-linear regime depend sensitively on the ini- 
tial perturbation. Thus setups Ml and M2 are used the follow 
the formation and evolution of a specific number of fragments. 

The time evolution of the equatorial plane density for setup 
Ml is shown in Fig|2] While the initial model is axisymmetric, 
it has already developed a strong m = 1 type deviation from 
axisymmetry at t = 6.43<d, which consequently evolves into 
a collapsing off-center fragment. At t — 7.45£d, we find an 
apparent horizon, using the numerical code described in 1 42] . 
The horizon is centered on the collapsing fragment at a coor- 
dinate radius of tah ~ 0.16r c , and has an irreducible mass 
of Mah ~ 0.24M sta r- Its coordinate representation is signifi- 
cantly deformed: its shape is close to ellipsoidal, with an axes 
ratio of ~ 2 : 1.1 : 1. The apparent horizon is covered by 
three refinement levels and 50 to 100 grid points along each 
axis. 

The evolution of the setup M2 is shown in Fig [5] Two or- 
biting and collapsing fragments are forming. However, even 
with the adaptive mesh refinement method we use, constraint 
violations prevent us from continuing the simulation to the 




FIG. 2: Time evolution of the equatorial plane density for setup Ml. 
Shown are isocontours of the logarithm of the rest-mass density. The 
four snapshots extend to 0.37r e and are taken at t/to = 0, 6.43, 
7.14, and 7.45, respectively. They show the formation and collapse 
of the fragment produced by the m = 1 instability. The last slice 
contains an apparent horizon demarked by the thick white line. Note 
that the shades of grey used for illustration are adapted to the current 
maximal density at each time, and that darker shades denote higher 
densities. 




FIG. 3: Time evolution of the equatorial plane density for setup M2. 
The snapshot times are the same as in Fig [2] In this case, two frag- 
ments are forming. Constraint violations have forced us to termi- 
nate the simulation before apparent horizons could be located, as ex- 
plained in the text. 
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formation of apparent horizons. Cell-based adaptive mesh re- 
finement, a better choice of gauge, or methods based on nu- 
merical analysis (e.g. I43I0 might be required in this case. 

To summarize, we have studied fragmentation and black 
hole formation of a general relativistic equilibrium polytrope 
of index N = 3. The polytrope has been shown to be unsta- 
ble to cylindrical perturbations of the form rsin(m</>), with 
rn = 1,2, which consequently grow to one or more self- 
gravitating fragments. We have applied an adaptive mesh 
refinement method to resolve the system accurately. In the 
m = 1 case, we have found an apparent horizon in the space- 
time, indicating that a black hole has formed. 

The dynamics of a nonaxisymmetric single-star collapse of 
this type is significantly different from that of the quasi-radial 
cases usually investigated. From the often considered case 
of quasi-radial collapse, to bar formation and subsequent col- 
lapse, to fragmentation and fragment inspiral we have a range 
of possible dynamical scenarios, which may be connected to 
discernable observable features in their gravitational wave sig- 
nature. In that sense, the evolution presented here can be con- 
sidered as an example of such processes. 

Extensive discussion of this work as well as further results 
for the binary case will be presented elsewhere. In this case, 
the fragmentation process could transform a star into a bi- 
nary black hole merger system with a massive accretion disk 
around it. 
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